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ABSTRACT 

The Navier-Stokes equations can be viewed as an incompletely elliptic 
perturbation of the Euler equations* By using the entropy function for the 
Euler equations as a measure of ^energy' for the Navier-Stokes equations, we 
are able to obtain nonlinear 'energy^ estimates for the mixed initial boundary 
value problem. These estimates are used to derive boundary conditions which 
guarantee L boundedness even when the Reynolds number tends to infinity* 
Finally, we propose a new difference scheme for modelling the Navier-Stokes 
equations in multidimensions for which we are able to obtain discrete energy 
estimates exactly analogous to those we obtained for the differential 
equation* 
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INTRODUCTION 


For computational problems involving the Navier-Stokes equations, it is 
necessary to limit the domain of computation and introduce artificial boundary 
conditions. Naturally, we would like these boundary conditions to be stable, 
compatible with weak boundary layers, and to remain valid even when the 
Reynolds number tends to infinity. Such a set of boundary conditions were 
proposed by Gustaffson and Sundstrom in [4], They used energy estimates on 
the linearized Navier-Stokes equations to obtain boundary conditions of 
maximal dissipative type. In this report we define an 'energy" in terms of 
the entropy function for the Euler equations and obtain fully nonlinear 
'energy' estimates from which we are able to extract a family of boundary 
conditions with the above properties. An attractive feature of these boundary 
conditions is that they are easy to implement and can be expressed in terms of 
the physics of the problem. 

The Navier-Stokes equations are an incompletely elliptic perturbation of 
the Euler equations — which are themselves a hyperbolic system of conserva- 
tion laws with entropy functions. It was observed by Mock [5] that by 
introducing the gradient of the entropy as a new variable a system of 
hyperbolic conservation laws can be reduced to a symmetric, hyperbolic system 
in terms of this new variable. Further, Harten [5] showed that if the 
dissipative terms in the Navier-Stokes equations are rewritten in terms of 
this new variable then the matrix coefficients of the dissipative terms have 
certain symmetry properties. We are able to show that the augmented matrix 
formed from these matrix coefficients is, in fact, negative semidef inite. 
This observation is crucial to the energy estimates we obtain for the Navier- 


Stokes equations 
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This leads us to propose a new difference scheme for modelling Navier- 
Stokes equations in multidiraensions . We are able to obtain discrete ^energy' 
estimates — which are exact analogs of the ^energy^ estimates we obtained for 
the differential equation — at the semidiscrete level, even for meshes with 
unequal mesh widths • Thus we are able to propose boundary conditions and a 
difference scheme for the Navier-Stokes equations which give a priori 
boundedness of "'energy' for all time. 

This report is organized as follows: In Section 2 we define the Navier- 

Stokes equations and obtain the necessary results to derive the 'energy' 
estimates of Section 3. In Section 4 we propose a family of 'stable' boundary 
conditions and relate them to the physics of the problem. In Section 5 we 
propose a new method for differencing the Navier-Stokes equations in 
multidimensions and obtain discrete 'energy' estimates for our difference 
scheme. Finally in Section 6 we obtain stable boundary conditions for the 
difference scheme and conclude by displaying some numerical simulations in 
Section 7. 

2. PRELIMINARIES 

We consider systems of hyperbolic conservation laws of the form: 

d 

(2.1) q^ + I f^(q)^ = 0. 

i=l i 

Here q(x,t) is an n column vector of unknowns, f^(q) is a vector valued 
function of n components, x = (x^,***,x^) and f = (f\**«,f^). 
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We can rewrite (2,1) in matrix form: 


( 2 . 2 ) 

where A^(q) 


i=l i 

The system (2.1) is called hyperbolic if the matrix 


(2.3) 


d i 

I w. A (q) 
i=l 


has real eigenvalues and a complete set of eigenvectors for all real 03 ^. 

Following Mock, a scalar function V(q) is an entropy function for (2.1) 
if; 

i) V satisfies 


(2.4) 


V f^ = 

q q q 


where F^(q) is some scalar function called entropy 
direction. 

ii) V is a convex function of q. 

T 

It follows from (2.4) upon multiplying (2.1) by 
solution of (2.1) also satisfies: 


flux in the x, 
1 


that every smooth 


(2.5) 




0 


where F = (f\... .F^^). 





The Euler Equations of Gas Dynamics 
Description of variables: 
p denotes density, 

u denotes velocity in the x direction, 

V denotes velocity in the y direction, 

w denotes velocity in the z direction, 

m^, and are the components of momentum in the x, y and z 

directions respectively, 

T is the temperature , 

p is the pressure, 

U is the thermodynamic entropy, 

E is the energy, 

R is the universal gas constant, 

Y is the ratio of specific heats* 

Note that we shall use (x, y, z) and (x^, x^) interchangably to denote 

jk 

the spatial vector x. 

We shall also need the following thermodynamic relations: 


T 


R 


P 


(y - 1) 




r 2 ^ 2 , 2^ 

m + m + m 
u V w^ 

2p 


U = log 



2^1 
m J 


- Y 


log p = log p - Y log P* 


up to an additive constant. 
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T, p and p will always be restricted to be positive because of obvious 
physical considerations, q will always denote the vector: 


p 


P " 

m 


pu 

u 


m 


pv 

V 


m 

w 


pw 

_ E _ 


_E _ 


The Euler equations are of the form: 


12 3 

q + f + f + f = 0 
^t X y z 


where 


■pu 


■ pv 


“ “ 




pw 

2 . 

pu + p 


pvu 


pwu 


^2 

2 ^ 

^3 


puv 

, f = 

pv + p 

, f = 

pwv 





2 

puw 


pvw 


pw + p 

(E + p)u 


(E + p)v 


(E + p)w 


We shall write the Euler equations in operator notation as: 


( 2 . 6 ) 


Eq = 0. 


The Euler equations have a family of strictly convex entropy functions defined 
by 

V(q) = -ph(U). 


The preferred entropy function in most physical applications is: 
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/ (m^ + + ra^) \ 

V(q) = -pU = -p log (e ^ —j + YP log p. 


The entropy flux functions turn out to be: 


/ r 2 ^ 2 ^ 2^ 

, / m + m + m 

f‘ = -m U = -m log E-— a ^ 

u u ° \ 2p 


+ ym^ log p 


F = -m U 

V 


F = -m U. 
w 


It should be noted that the entropy function V(q) is strictly convex but may 
be nonpositive in general. 


Navier-Stokes Equations 

We shall denote the Navier-Stokes equations in operator notation as: 


(2.7) Nq = 0 

where 

Nq = Eq + (-D)q, 


where (-D)q 



3 

I A^J(q)q 
j = l j 


X. 

1 


We can represent (2.7) in the alternative form: 


3 

Nq = Eq + I (h^)x. 

i=l ^ 



where 


3 

h"= I 

j=l 


A^^(q)q 


X. 


12 3 

and h , h , h are as follows: 




0 
e 
e 
0 ^ 

ax 

9 u + 9 v + 0 w-k — 
XX yz zx 3x 


XX 


yx 


= 


0 
0 

} 

S 

3x 

0 u+9 v+0 w-k— 

xy yy zy 3y 


xy 

’yy 


h’'. 


xz 


yz 


zz 


0 U + 0 V + 0 W-k — 
xz yz zz 3z 


where: 


®xx - ^ ~ ^ ^ 
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Here the parameters y, X and k are as follows: 
y is the shear coefficient of viscosity, 

X is the second coefficient of viscosity, 
k is the coefficient of thermal conductivity. 

Clearly 

12 3 

(-D)q = h + h + h . 

X y z 

We write the Navier^Stokes equations as 
(2.8) Nq = Eq + (-D)q 

where E is the differential operator corresponding to the Euler equations 
and D is the elliptic perturbation (depending on the parameters y, X and 
k) due to the dissipative terms in the Navier-Stokes equations. 

The entropy function V(q) = ~pU for the Euler equations is not positive 
valued in general. Since we wish to interpret the entropy as a measure of 
'energy' of the system, we can normalize V(q) and define a new entropy V(q) 
which has the properties: 

(i) V(q) > 0 Vq 

(ii) V(q) = 0 <=> q = q for some fixed q. 

For the entropy V is not altered by adding to it an arbitrary inhomogeneous 
linear function. Hence we define V(q) as follows: 
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(2.9) V(q) = V(q) - V(q) " I (q)(q. " q.). 

i=l ^ ^ 

The associated entropy flux functions are given by 

(2.10) F^(q) = F^(q) - F^(q) - I (q)(f^(q) - f^(q)).. 

We choose as our rest state q = (o, ra , m » m , E) such that: 

^ ^ u V w 

7 > 0 

u = V = w = 0 
T > 0. 

In a later section we shall put further restrictions on p and T, With this 
choice of q we obtain: 

V(q) = p U - U + + w^) - Y Z ■*" (y “ Op 

[ 2RT T J 

N _ “ TT . (y "" i) / 2 , 2 , 2,, , T 

F (q) = pu U - U + — i — — — (u + v + w)-Y"*"” • 

L 2RT T 

Or more compactly: 

F^(q) = u(V(q) - (y - D'p) 

F^(q) = v(V(q) - (y - l)p) 

F^(q) = w(V(q) - (y - l)p). 
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We make the change of variables v 
terms of v. Here 


and rewrite the operator Dq in 


V = 


a 

b 

c 

d 

e 


Substituting these relations in 


(-D), = j ( j \ 

1=1 \ J=1 j/ 1 


yields 


3/3 


(-D)q = (-D)v = W I C ‘^(v) 
i=i \j=i 


X. / X. 

J/ 1 


where C^-^(v) = A^"^(q)q • Harten [5] observed that the matrix coefficients 
C^^(v) satisfy the symmetry relation: 


For the energy estimates we shall derive in Section 3 we need to show that 


( 2 . 11 ) 


3 3... 

I I (5^)^ C^J(v)(5h < 0 ¥ e I^. 

i=l j=l 


Clearly this is equivalent to proving that the augmented matrix C defined by 
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r^ii oi2 

c c 

a= 

is negative semidef inite. 

Observe that C is symmetric. Hence to check that C is 

A 

negative seraidefinite it is enough to perform an LU decomposition of C and 
then check that the diagonal elements of the upper triangular matrix U are 
negative. Thus we obtain: 

(2.12) (2.11) holds <=> X + li >_ 0. 

Typically, for physical fluids ^ ^ ^ y ^ 0. Further for most 

fluids under nonextreme flow conditions the following relation holds: 

^ X + 2p y. 

Here Pr denotes the Prandtl number. So (2.11) holds under these conditions. 
3. NONLINEAR ENERGY ESTIMATES 

Let ^ denote a bounded domain in 1^ and let 9^^ denote the boundary 
of fi. Consider the mixed initial boundary value problem in Q: 

(3.1) Nq=0 VxeO, t>0 



where N is the Navier-Stokes differential operator with initial condition: 
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(3.2) q(l^, ^ e U 

and boundary condition: 

(3.3) Bq = 0 V X e an, t > 0 


where B is a boundary operator. 
Define: 


(3.4) S(t) = / V(q(x, t))d^. 

n 

We claim that S(t) gives us an estimate of the energy of the system at 
time t. Note that V(q) is a strictly convex, non-negative function of q, 
i.e. , 

i) V(q) > 0 

ii) V(q) = 0 <=> q = q 

iii) Vqq > cl where c > 0, at least in some appropriate physical domain. 
In particular, (iii) implies : 

(3.5) V(q) >_ allq - qll^ 

for some a > 0. Hence we conclude that S(t) has the following properties: 

i) S(t) > 0 

ii) S(t) = 0 <=> q(x, t) = q V x e Q 

iii) S(t) < k, where k is a constant => 
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|/ nq(x, t) - qll^ h 


if q(x, t) lies in a domain where an inequality of the form (3.5) holds. 

We wish to study the time evolution of S(t). Recollect that the Navier- 
Stokes equations are: 

Nq = Eq + (-D)q = 0 


where 


3 3 

(-D)q = I I 
i=l j=l 


J 


i 


Making the change of variables v = V 

q 


3 

(-D)q = (-D)v = \ 

i=l 


gives : 


1, 

J = 1 J 1 


where the matrix coefficients C^J(v) satisfy the condition: 


(3.6) 


I I C^J(u)5J < 0 ¥ 5^, € ]^ . 

i=l j=l 


This leads us to the following theorem. 


Theorem 3.1: Consider the mixed initial boundary value problem (3.1)- 

(3.3). Let denote the outward unit normal to the boundary 9Q 

and let f\ • • • denote the entropy flux functions as in (2«10). Then any 
piecewise smooth solution of (3.1)-(3.3) satisfies the energy estimate : 


(3.7) 


f </ I 


I I 


3 3 ... 

V A^J(q)q 

i=l j=l ^ 


da. 
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Here da is an element of surface area. 


Proof ; We prove it here for the case q(x, t) is 'smooth enough.' We 


have 


(3.8) 


Nq = q. + I f^(q)„ + I ( I A^'^(q)q ) = 0 . 

i=l ’'i i=l \ j=l ""j / ""i 


Premultiplying (3.8) by we get 


\ j, A^j(q)q^_ I = 0. 


i=l 




By (2.4), this reduces to 


(3.9) 


3 . 3 / 3 . . 

\ + I + I V ( ^ A^J(q)q 

1=1 1=1 ^\j=l 


^j/^i 


= 0 . 


From (3.9) we obtain 


dS 


1=1 ^i i=i Mj=i j/i 



dx. 


or 


f--/ ''M j, 

Q ( 1=1 1 1=1 \ J = 1 


"j/"i 


dx. 


By the divergence theorem 


dS 

dt 


/ 3 . . 3 3^. 

- / I I I 

dSl \ i=l i=l j=l 


(v)v 1 da 

^j/ 



dx 


+ / 


3 3 


I I 

1=1 j=l 


v'^ C^^(v)v 'j 

X. X. / 


But by (3.6) 


(A \ <■ 


dx < 0. 


Hence we obtain 


dS 

dt 


< 



3 . . 3 3 

I I I 


i=l 


1=1 j=l 


v'^ C^^(v)v^ 

j 


da* 


This finally yields 


(3.7) 


jc / 3 . . 33 

§< - / I f I I 

\i=l i=l j=l 


V A^'^(q)q 
q ^ ^ ^x. 

J 


da. 



Remark; Estimate (3.7) is a fully nonlinear ''energy' estimate which 
holds for a broad class of solutions of the mixed initial boundary value 
problem (3.1)-(3.3). In fact (3.7) is simpler to obtain and broader in scope 
than the linearized energy estimates in current use. 
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4* STABLE BOUNDARY CONDITIONS FOR THE NAVIER-STOKES EQUATIONS 



Figure 4.1 


Since the Navier-Stokes equations are rotationally invariant, we can 
consider a moving coordinate frame (x, y, z) where x points in the 

direction of the inward normal and y and z are tangential to Of 

course, we reorient u, v, and w so that u is the component of velocity in 
the X direction, v in the y direction, and w in the z direction. Let 
^ denote the outward unit normal to 3Q. Clearly ^ = (-1, 0, 0). 

Then the nonlinear energy estimate (3.7) takes the simpler form: 


(4.1) 


\ J=1 


V q 

q X , 

J 


dc 


or 
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(4.2) 


^ < / (f^ + V hhdc 

St - '• q ^ 




where 


h‘ - I » . 


. , "X. 

J=1 J 


Substituting the relations obtained in Section 2 we obtain 



So finally we obtain the needed 'energy'' estimate: 
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(4.4) 


dS , f , f k(Y - 1) rT - 3T ^ 

— < / u(V - (y - l)p)dCT + / ^ I— r— ) ^ da 

3n 3J1 RT T 


-/ v„(|i . |i)d, - / * |i)d< 

3t2 RT ^ 3fi RT ^ 


where 


(4.5) u(V - (y ~ Dp) = pu 


- „ . (y - 1) / 2 ^ 2 , 2^ , T 

U - U + (u + v + w)-Y+ — 

2RT T 


and 


(4.6) V = p 


U - U + h-l-ll (u^ + + w^) - Y + - 


2RT 


+ (y “ Dp. 


So far we have left T and p vaguely defined. We are now going to 
specify T and p , but first we have to make a few assumptions about the 
solution q(x, t) to the mixed initial boundary value problem (3.1)-(3.3). 

Assumptions: 


i) 3 p . > 0 3p(x, t)>p. VxcQ, t>0 

^ ^min K > y _ '^rain ’ — 


(4.7) ii) 3 p > P • ^ p(x, t)<p VxgQ, t>0 
^ ^max ^rain ^ ' ^max * — 


iii) 3 >03 T(x, t) > V X e t > 0. / 
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Assutnptions (i)-(iii) exclude cavitation, freezing of the fluid, and such 
esoteric phenomena as the formation of black holes. Further, it should be 
noted that if (4,7) is valid V(q) is strictly convex and hence, then we are 
able to obtain boundary conditions which ensure that S(t) remains bounded, 
i.e, , 

S(t) ^ k V t _> 0 where k is a constant; then by the discussion 
following (3,5) 

Hq(x, t) - qll^* die C (another constant) V t _> 0, 

2 

Hence we obtain L boundedness of the solution for all time. 

Choose 0 < p < p , such that 
min 

(4.8) V - (y - 1)7 > 0 


and 0 < T < T . such that 
rain 


(4.9) 


U - U < 0 


V q staisfying (4.7). 

It is easy to see that this can always be done. For by (4.6) 


or 


V(q) > p 


U - U - Y + Z 

T 


+ (y " l)p 


V(q) > p 


(y - 1) log (-2.) ~ Y + z “ t) 
P T T 


+ (y ” l)p- 
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Clearly 

^ - log (z) > 1 ^ T > 0. 
T T 


Hence we get 


V(q) > p 


(y - 1) log (z) “ Y + 
P 


+ (y “ i)p» 


Thus if p > p . >0 we can choose 


p < p . such that 


V(q) - (y “ l)p > 0. 


Hence (4.8) holds. 


Now 

U-U = (y- 1) log (^) - log . 

. P T 


If (4.7) is valid then 


U - U < (y - 1) log (-^5^) - log 

P j 


We can choose 0 < T < T . 

min 


(y - 1) 

P 


such that 

P T . 

. ^^max^ _ r rainN 
log [— zr-J “ log < 0* 


Hence (4.9) holds. 

The Navier-Stokes equations are: 


(4.10) 


Nq = Eq + (-D)q = 0 



- 21 - 


where E is the hyperbolic differential operator corresponding to the Euler 
equations and D is the incompletely elliptic perturbation, depending on the 
parameters, y, X, and k, due to the dissipative terms in the Navier-Stokes 
equations. It is assumed that y, X, and k are proportional to a small 
parameter e > 0. Thus the Navier-Stokes equations may be viewed as an 
incompletely elliptic perturbation of a system of hyperbolic equations. 

The question we are concerned with is: for which boundary conditions are 
the solutions of the perturbed problem (i.e., e > 0) well defined in some 
time interval 0 ^ t 'Tq and are bounded in an appropriate norm as e ^ 0. 

Michelson [8] suggests boundary conditions of the type given below for 
the singular perturbation problem: 

(4.11a) Bq = S( { (eD)“}^^^; x edil, t > 0, e > O) = 0 
where 



and a = (a^ a 2 > a^) is a multi-index in some finite set A _c it . 

It is well known that a singular perturbation problem of the type we are 
considering exhibits a boundary layer phenomenon. We wish to choose boundary 
conditions in such a way that a strong boundary layer does not develop. For 
the boundary conditions to be compatible with a weak boundary layer, we need a 
condition of the sort given by: 

®le=0 t 0, e = O') = 0 


(4.11b) 
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gives a well-posed boundary value problem for the hyperbolic part of the 
Navier-Stokes equation: 

N|,.0 I - E, - 0 

(see [8]). 

Further, from Strikwerda's work on initial boundary value problems for 
incompletely parabolic systems, we know that the number of boundary conditions 
which should be imposed for the Navier-Stokes equations to obtain a well-posed 
problem are: 

5 boundary conditions for inflow boundary and 

4 boundary conditions for outflow boundary. 

At the same time the unperturbed hyperbolic system requires: 

5 boundary conditions for supersonic inflow, 

4 boundary conditions for subsonic inflow, 

1 boundary condition for subsonic outflow, and 
none for supersonic outflow. 

The boundary conditions we are going to impose will be of the form: 

(4.12) e R |a + Sq = g 

where R is a matrix of rank at most 4. 

To get a set of boundary conditions that also works for e = 0 S must 
be chosen in such a way that Sq = g is a proper set of boundary conditions 
for the unperturbed hyperbolic problem. If Sq = g gives too many boundary 
conditions for the unperturbed hyperbolic problem, the solution will contain 
boundary layer components of the form exp(-x/e) for e 0; (see [4]). 
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We now proceed to suggest a set of maximal dissipative' boundary 
conditions for the Navier-Stokes equations which are compatible with weak 
boundary layers. 


Outflow Boundary Conditions 

We need to specify 4 boundary conditions for the perturbed problem. For 
supersonic outflow (-u > c > 0) we need not specify any boundary conditions 
for the unperturbed hyperbolic problem, while for subsonic outflow (c > -u > 
0) we need to specify only one. From (4.4) we have 


y ( f*\i f i\ — \ j - I r ^(y 1) /-T "" 3T 

— < u(V - (y - l)p)da + -3; J ^ da 


3a 


3a RT 


3x 




t (y “ 1) f (y - 1) r9w . 


By (4.8) 


f u(V - (y “ l)p)da £ 0. 

3a 


The boundary conditions we suggest which give a decay of energy are: 
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(4.13) 


k(Y - 1) 3T . „ T = ^ 
R 3x 1 ^ 




f^v ^ 3u-> 

37^ 


“3 ^ 


= g. 


/3w , 3 un 
'^Ux 3z^ 


«4 w 


S/ 


“l 2. ^ 

«2 u = g2 , a2 > 0 

03 > 0 

“4 > 0 * 


It should be noted that if we put = 0 we must also put = 0 for 

i = !,•••, 4. Since we want (4.13) to be compatible with the preceding 

discussion, we end up with the following set of maximal dissipative' boundary 
conditions: 


Supersonic outflow 


k(Y - 1) ax ^ Q 

R 8x 




= 0 


(4.14) 


rdu 


+ 


9v 

3x 


) - » 



+ - 0 

dz^ 
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Subsonic Outflow 

k(Y - 1) 9T ^ Q 

R 8 x 


( 4 . 15 ) 


C2p + X)f + 



u = gn 


02 > 0 


,(|v + |H) = 0 

^''9x ay-* 


f^w , 9u., 


Inflow Boundary Conditions 
By ( 4 . 9 ) 


/ pu[U - U]da £ 0 
9 fl 


Subsonic Inflow 


Since 


c > u we have 


2 

u 


< yRT. 


Hence 


J (pu)u^ S. ! ~ — (pu)YRTdo 

9 fi 2 RT an 2 RT 


. dS . , (y - 1 ) //o . ,N 9 u , , aw-.,, 

=>dFl"/ u ( 2 y + X) - 5 ;^+ X(— + ^llda 

SO RT V 


3 x ^ 3 y 3 z 


r (y “ 1 ) / (pu) ^ 

+ / — rr— ''^—r " - VV'*” 


an RT 




+ / = 
an T 


( I - J.) - .1) II + iT . -- - T I 2 ) 

T R ax 2 ^ 


da. 
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The boundary conditions we could specify are; 


pu = g, 


If * '(I7 " ■ “ 


(4.16) 


fdv . 9u-> 

" “3 ^ = §3 


rdW dU\ 

‘'fe " “4 ^ “ §4 


MXJLilH- „ T = g , 

R 9x “5 ®5' 


Here we must impose the conditions 


“3 ^ 

«4 > Si/2 

(y^ “ y 4” 2) 

> i ^ ^ some suitable e > 0 


depending on T 
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Supersonlc Inflow 


dt 


< / pu[U - U]da + / u 

3n 9fi RT 


(pu) 

2 




• 9v 9w ^ 

'9y 9 z- 


da 


, . (y - 1) 

+ J - - V 

9n RT 


(pu) ^ 9 v , 9u->j_ 

~r ■ "fe * 


+ I 


(y - 1) 


w 


RT 




da 


+ / = 

9n T 


^ ^ i ^ 


da. 


The boundary conditions we could specify are: 


(4.17) 



g 


2 


k(y - 1 ) ^ 
R ax 


05 T . 


Here the following conditions must hold: 


a 


2 ’ 


“4 



and + e for some suitable e > 0 
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Dlscussion of the Boundary Conditions 

The outflow boundary conditions we have proposed are of the following 

type: 

(4.18) -F = 0- 

9x 

(4. 18) says that the computational boundary corresponds to an insulated wall 
and there is no conduction of heat across it. 


(4.19) 


XX 


- 0 


yx 


zx 


0 


(4.19) asserts that there is no shearing of the fluid either in the normal or 
tangential direction against the computational boundary. 

For an inflow boundary one of the boundary conditions we have specified 
is: 


(4.20) 


pu = g. 


In other words , we specify the momentum influx across the computational 
boundary. 

For the temperature component of the state vector, we have the boundary 
condition: 


k(Y - 1) ^ _ 


9x 


«5 T = g5 


(4.21) 
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(4,21) has an interesting physical interpretation. For injection through a 
porous wall into the main stream (sometimes called transpiration) the injected 
fluid may be, say, a coolant at a temperature considerably different from the 
wall, and one needs to consider an energy balance at the wall, A good 
approximation for coolant injection is to use the boundary condition proposed 
by Roberts [in Truit (1960, chapter 11)]. 


(4,22) Injection: 


k|-=puc(T-T . J 
3x ^w w p w coolant 


where u^ is the momentum flow of coolant per unit area through the wall, 

is the temperature of the wall, T^oolant temperature of the 

3T 

coolant and — ^ is the temperture gradient at the wall. (For more 
information refer to [12, chapter 1,4],) 

By choosing and g 5 appropriately (4,21) can be geared to satisfy 

(4,22). The physical effect of such a boundary condition would be to make the 
temperature of the fluid within the computational boundary stabilize at 
Tcooiant time. When we discretize the boundary conditions, however, we 

can choose so that the discretized version of (4,21) is of extrapolation 

type. This has the effect of allowing the system to evolve as if there were 
no boundary. 

For the velocity component of the state vector, the boundary conditions 


we have specified are: 
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(4.23) 


(2y 


+ X) — + Xf— + — ) 


'9y 9z’ 


a 2 u = §2 



V = gq 


w = 


§4* 


For one-dimensional fluid flow the boundary conditions (4,23) would take 
the simple form: 

(4.24) (2y “ “2 “ ®2* 


For perfectly diffuse reflection of a gas or fluid against a solid wall 
we have the boundary condition (after Maxwell in 1879): 


(4.25) 


w 




8x''w 


where is the velocity of the fluid at the wall, £ is the mean free path 
of the gas or fluid at the wall and (|^)^ velocity gradient at the 
wall (see [12]). 

By choosing g 2 and (%2 appropriately we can choose (4.24) to satisfy 
(4.25) . Since we want the boundary conditions to be nonreflecting, however, 
this would be a bad choice indeed. If we let the computational boundary move 
with the same velocity as that of the fluid at the boundary interface, we 
should not have any diffuse reflection. And, in fact, we can discretize 
(4.24) to achieve this exactly. The discrete boundary conditions so obtained 
are of extrapolation type and the effect of the computational boundary is 


minimized. 
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A last remark is that the boundary conditions proposed are of maximal 
dissipative type and such boundary conditions are intrinsically radiative* 


5. STABLE SEMI-DISCRETE DIFFERENCE SCHEMES FOR THE DIFFERENTIAL EQUATION 

The stable difference scheme we are going to propose is valid for the 
Navier-Stokes equations in multidimensions. However, to avoid getting drowned 
in notation we restrict the discussion that follows to Navier-Stokes equations 
in one space dimension. 

Consider the incompletely parabolic system obtained from a hyperbolic 
system of conservation laws with entropy: 

(5.1) q^. + f(q)^ + 

Let V(q) be a normalized entropy function as in (2.9) for the 
hyperbolic system of conservation laws: 

(5.2) q^ + f(q) = 0 

t X 

T 

and let v = V . 

q 

We can rewrite (5.1) as: 

(5.3) q + f(q) + (C(v)v ) = 0 


where we assume C(v) < 0. 
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Divlde the x axis into a mesh of points {x.} where x. = iAx 

and Ax is the mesh width* 

Let Consider the seraidiscrete difference 

approximation to the hyperbolic system of conservation laws (5.2); 


(5.4) 


dt 


A" 


— (h(q. ,***,q. . ) - 0 

Ax j-r* j+s 


where r and s are integers >_ 0 and h(q^_^, • • • ,qj_j^g) is the numerical 
flux function corresponding to the flux function f(q), i.e., h(q. >***>q.. ) 
satisfies the consistency condition: 


(5.5) 


h(q,***,q) = f(q). 


Let the order of accuracy of the semidiscrete difference approximation be a 
where we can take a >_ 1. 

Further, suppose (5.4) satisfies a semidiscrete entropy inequality: 


(5.6) 


dV(q^ ) 

dt 


[H(q._^,---,qj+3)] < 0 


where H(q._ >***>Q-. ) is the numerical entropy flux function corresponding 
j— r J'^’S 

to the entropy flux function F(q), i.e., H(q._ j***>q.. ) satisfies the 
consistency relation: 


(5.7) 


H(q, • • • ,q) = F(q) . 
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We propose the following semidiscrete approximation for the Navier-Stokes 
equation (5.1): 


(5.8) 




’Vs” 


1/2 ’ 


A.v. 

7 ^) = 0 

Ax ^ 


where 

V. + 

1/2 ■ ^ 

It is easy to see that the order of accuracy 3 of (5.8) is given by 
3 = rain(a, 2). Hence if (5.4) is a second-order accurate approximation to 
(5.2) then (5.8) is a second-order accurate approximation to (5.1). In any 
case, (5.7) is at least a first-order accurate approximation to (5.1). 

Remark ; The only second-order accurate semidiscrete difference scheme 
for a hyperbolic system of conservation laws which satisfies an entropy 
condition of the type (5.6) of which we are aware of is one proposed by Osher 
in [10]. 



Lemma 5.1: The semidiscrete difference approximation (5.8) to the 

Navier-Stokes equations (5.1) satisfies the entropy inequality ; 


vT 


(5.9) ^ (q.) + ^ [H(q._^,...,q.^g)] + A - (C(v < 


0 . 


Fix At > 0. Let 


At 


qMt) - q^ - _ A - (h(q^_^,...,q^^^)). 


Let q (t) denote the solution to the semidiscrete difference scheme 
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(5.10) 



^ A- . . E E ^ 

Ax ^^j-r * * * * ’ '^j+s ^ 


= 0 


F ^ F 9 

with ^ j • Then q^CAt) = q^CAt) + 0(At ). Since (5.10) satisfies the 

entropy inequality: 




V sequence 0 3 a sequence {e^} 9 ->• 0 and the discrete 

entropy inequality: 


(5.11) 


At, Ax 




i 'k 


holds. 

We now define a discrete difference approximation for the Navier-Stokes 
equation: 

q.(At) = q*(At) - A - (C(v . 


Let denote the solution to the seraidiscrete difference scheme: 


dq . 

^ A- 


(h(q ,...,q )) + — fC(v ,, ) —^—^1 = 0 

dt Ax ^j“r Ax ^ j+ V 2 Ax ^ 


A. V. 


Then it is obvious that q^(At) = qj(At) + 0(At ). Define an approximate 
entropy function 


T At 


^ -j. 


W(q.(At)) = V(q.(At)) - v . — A - (C(v.^ ) — ^) 
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By (5.11) 


At, 


W(,j(At^)) < V(,j) - _4 - 


- vj ^ (C(Vj, V2 > “k- 


or 


(5.12) 


U(? (it ) - V(q ) 

■ - - it : — - If 


T 

V. 

_JL 

Ax 


1/2 ^ 


By a Taylor series argument: 


W(q^(At)) = V(q^(At)) + O(At^), 


And since q^CAt) = qj(At) + 0(At ) => 


T 

V , 


^ <ij> < -Ii t«Vt--'V3>l - ii‘ - 


by letting ^ 0 in (5.12). 

Define a discrete version of the energy S(t) by: 


N 


(5.13) 


S(t) = I AxV(q.), 

j = l ^ 


We can now easily obtain a discrete version of the energy estimates of Section 
2. (Henceforth we shall denote ) by hj and 

H(qj_^.---,qj+ 3 ) by H^.) 
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Theorem 5,1: Consider the semidiscrete difference approximation: 


(5.8) 




dt 






and let 


S(t) = ^ AxV(q ), 

j=l ^ 


Then the following ^energy^ estimate holds; 


(5.14) 


“O ^ ''l 1/2 ' — 


^ ''n 


^ ''' ^N+1 V 2 ^ Ax 


Proof: By Lemma 5,1 


dV(q.) . ™ A, V. 

< - ST IHj j - V. - [C(Vj^ 1/^ ) -jji] 


Hence 


N dV(q.) N 


f .1 n - iHj 1 - .J, ''j 1/2 > 

j=i j=i j=i j ^ 


=> ii 


4 + Vj C(Vj^V2> 


^ '' .1 
Ax 


T 


A V 
+ N 


^1 1/2 ^ ~sr ■ \+i 1/2 ^ ~ 


using summation by parts. Since 


N A^ V 

C < 0 -> I i + V. ) — ^ < 0 
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Hence 


< + V, C(vi,_ ) 




T 




Remark ; The generalization of this theorem to multidimensions is 
obvious. 

We now extend these results to uneven meshes, though restricting (5.4) to 
the Godunov scheme. 


The Godunov Scheme 


Consider a grid of points and define Ax^ 


- Xy Define 


X. 1 + X. Ax. , + Ax, 

X. 1 , 7 z and ox. r 

J~ V2 2 J 2 



Figure 5.1 
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Conslder the hyperbolic system of conservation laws: 


+ f(q)x = 0 


(5.15) \ with Rieraann initial data; 


q(x, = 


q. for X < 0 
J 


*j+l for X > 0 


G G 

Then the Godunov flux f (q . , q . , , ) = f . is defined as; 

J ^J+1 J 


(5.16) 


f^(qj> qj+p = f(q(j^ = o, t = o'*’)) 


where q(x, t) is the solution to the Rieraann problem (5.15). 
The seraidiscrete version of the Godunov scheme is: 


(5.17) 


A- r 


where f (^j > ^j+1^ defined in (5.16). The seraidiscrete Godunov scheme 

satisfies the entropy inequality: 


(5.18) 






G G 

where F (Qj > numerical entropy flux for the Godunov scheme 


F^(qj, qj+i) = F(q(x = 0, t = o"^)) 


and is defined as: 
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where q(x, t) is the solution to the Riemann problem (5.15). 
A natural modification of (5.8) for uneven meshes is: 


(5.19) 


^ ^ fij ^ ^ V 1/2 ' “ 


0 . 


Clearly (5.19) is a first-order accurate semidiscrete difference approximation 
to (5.1). We claim that (5.19) satisfies the entropy inequality; 


(5.20) 


A- , 

dt ^ ’ ‘^j+l^^ 


+ V 


,T ^ 

j 


(C(^ 


j+ V 2 


V 


Ax. 

J 


i) < 0. 


The proof is exactly the same as in Lemma 5.1 • 

Define a discrete version of the 'energy^ S(t) for the uneven mesh 
{Xj} by: 


N 

(5.21) S(t) = I 6x. V(q.). 

j=l ^ ^ 


Then we obtain the following theorem. 


Theorem 5.2: Consider the semidiscrete difference approximation : 


(5.19) 


dqj = _ ^ 

dt 6x 



kL 

6x. 


(C(v, 


j+ V2 


A + V. 


N 

S(t) = I 6x V(q ). 

j=l ^ ^ 


and let 
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Then the following ^energy^ estimate holds; 


(5.22) 4^ < 

at “ 


''l 1/2 ’ 


n -I- ''O 

Ax„ 


f:, + V, 


N+1 


« V 1/2 ’ 


Ax,, 


Proof ; The proof is exactly as in Theorem 5.1. 

Remark ; Theorem 5.2, though a simple extension of Theorem 5.1, will prove 
very useful in proposing stable boundary conditions for seraidiscrete 
difference approximations to the Navier-Stokes equations. 
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6, STABLE BOUNDARY CONDITIONS FOR THE DIFFERENCE SCHEMES 



Figure 6.1 


Suppose we have the boundary condition at x = 0: 


g T = g where 8 can be arbitrarily small. 

oX 


Then a discrete version of this boundary condition would be: 


(T - T ) 

Ax„ 


- Tq = g 


8T^ - gAxQ 
^0 " (e + AXq) • 


So if AXq » 3 => Tq “ g and hence could give rise to a steep ^numerical 
boundary layer.' To avoid this we should have AXq of the same order of 
magnitude as 3 ov smaller. If, however, we choose an even mesh this becomes 
computationally infeasible since 3 can be arbitrarily small. 

To overcome this difficulty we propose the following: Divide the 

interval [xq, into N + 2 points Xq, before let 

= *j+i - ’'j- 


Choose 
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Ax^ = Ax-_ = Ax^ 

0 N 

Ax. = Ax for i = N-1 

1 

where we may choose Ax^ « Ax. (See Figure 6.1.) So the mesh points are 

evenly spaced in the interior, but xq is close to x^ and is close 

to Xjj. 

Ax "t* Ax^ 

Define 6x ^ • Henceforth we shall restrict q to be 

T 

q = (p» E) • 

The semidiscrete difference approximation for the mixed initial boundary 
value problem is tailored according to (4.19): 


dqi 

+ ^ 

dt 

6x 


+ ^ 

dt 

Ax 


+ Az 

dt 

6x 


A + v^ A + Vq 

3/2^ ~A^E ^ Ax' 


C(v„,„) 


= 0 


A- 


A + V, 


V 1/2 ' ~ 


* Rv ^‘In'Ih+1 q Sir 


j 

A + V, 


= 0 for j=l , • • • ,N-1 


N 


A + V, 


1/2 * V 1/2 ^ JT 


N-1 


= 0 


( 6 . 1 ) 


with initial condition: 


and boundary conditions: 


qj(0) = Qj 


B(qQ, qj, q2> = 0 Vt 


“(%_1» %> <1n+i) 0 Vt 


where B and B are boundary operators. 
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In general, B and B may be an underdetermined set of boundary 
conditions. For example, b(<^Q> ^ 2 ) form: 




T = T 
0 1 


and p 


0 


is unspecified. 


If B is a ^legitimated boundary operator we can connect to qg by just 

two waves, say a 2- contact and a 3- shock or rarefaction, such that 
S(qQ, q^y q 2 ^ ~ ^ holds. We make this statement more precise by examining 
the boundary Rieraann problem. 

Suppose we have a boundary (a wall) given by x = st. We want to give a 
number, m, of nonlinear boundary conditions depending on q at the boundary: 




h: IP + iP, heC” 


( 6 . 2 ) 


and initial condition: 


q(x, t) I = q^ (a constant state) for x > 0. 


The underlying differential equation is 


( 6 . 3 ) 


A(q)q^ = 0 


where 
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A(q) 


af 

3q 


We want to construct similarity solutions for this initial boundary 
problem. Clearly, for the problem to be well-posed we have to choose q^ so 
that 


where Xi(q) are the eigenvalues of the matrix A(q) and m = n - k* 

We can connect q^ to the boundary by n - k waves. The raref unction 

th 

and shock curves of the j — family determine a wave curve: 

e.> 0 

Vs^(q^, t.), e. < 0 

where is twice continuously differentiable in both its arguments (see 

[ 2 ]). 

Given q^ and sufficiently small parameters ^k+l*****^n define 

states q^,»**,q^^^ inductively by: 


Set: 


q = W (q^, e^) 


qj ^ = W*^ (q*^ , e^. ) for j = n, •••jk + 2, 


_ '^k+1^ V 
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Then the boundary condition h(q) = 0 becomes ra equations in the n - k 

unknowns e, .,,•••»£ • 
k+1 ’ n 


By the implicit function theorem we will be able to solve for in 

k+1 ’ n 

r\ Vl 

terms of q- in a manner if ra = n - k and the Jacobian has rank 

^ dS 

m. Differentiating we get the local solvability condition: 


(6.4) 




has rank 


m 


n - k 


where r^ in the i— eigenvector of the matrix A(q) corresponding to the 
eigenvalue The construction is described in Figure (6.2). 



Figure 6.2 
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(6.4) is also a necessary condition for the linearized initial boundary value 
problem to be well-posed. 

We relate this now to the boundary operator 
denote and denote q^ , Then <^2 ^ ® corresponds to 

the ra boundary conditions: 

(6.5) ° 

and the initial condition for the initial boundary Rieraann problem (6.2) is: 

q(x, " 'Ir " 

So B(q^, q^, q 2 ) is a ^legitimated operator if (6.4) holds. Assuming that 
it is, we can connect q^^ to a state q^ such that ^2^ ~ ^ 

holds • 

Returning to the difference scheme (6.1) we can now define the Godunov 
flux f (qg, qj^) using this construction. Suppose, at t = 0 qj^(x, t) = q^ 
and we have the boundary operator ^2^ ~ 

Consider the Riemann initial boundary value problem: 

q(x. for X < 

and boundary condition: 

h(q) = 0 for x = Xq corresponding to B(qg, q^ , q 2 ) = 0 


as in (6.5). 
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Figure 6.3 

Then if (6.4) holds we can join to a state qq in a small enough 

neighborhood of so that BCq^, q^ , q 2 ) = 0. 

In general, the number of boundary conditions, m, specified by the 
boundary operator B will be more than the number of boundary conditions 
needed for the hyperbolic Euler equations since we are solving the Navier- 
Stokes equations, i.e., ra > n - k. Hence, when we connect to a state 

qQ lying to the left of the line x = x some of the waves will lie in the 

region x < x . This corresponds to letting waves be radiated at the 

boundary. 

We know that we can connect q^^ to a state qg such that > 

q^) = 0 if qQ and q^^ are sufficiently close. We chose an uneven mesh 
precisely for this reason. We set Ax^ = Ax^ = Ax^ and Ax^ = Ax for 
i = !,•••, N - 1 where Ax' could be arbitrarily small. In fact, by choosing 
Ax' small enough we can make qQ lie as close to q^ as we want. This has 

the effect of making the waves in the cell bounded by x = xq and x = X|^ 

very weak and hence we conclude that 


F^(qQ» ^i) = qp = FCq^)' 
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Recollect that the energy estimate we obtained for an uneven mesh was of 
the form: 

f <[fX- 


( 6 . 6 ) 


A + V, 


\+i I/2 ^ Ax' ] 


Since qg ->■ as Ax' -► 0 the estimate (6.6) reduces to: 


(6.7) f < [¥(q,) + V (,^)A(,j) - [F(,^) + V (,^)A(,^) 


for sufficiently small Ax'. 

One further point that should be noted is that when we take a fully 
discrete version of our semidiscrete scheme the CFL condition we have to 
impose to prevent wave interactions is: 


( 6 . 8 ) 


AtX(A) < 


Ax + Ax' 
2 


where X(A) is the spectral radius of the matrix A(q). Clearly (6.8) is not 
an unduly restrictive condition. 

For the numerical simulation presented in the next and last section, we 
nondimensionalize the Navier-Stokes equations. 

Let L be a reference length, a reference density, and u^ a 

reference velocity. Define : 

* 

X = x/L, 

* 

t = tu^/L, 

* 

u = u/u^, 
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* 2 

P = P/Pf u^, 

* 2 
T = TR/u^ , 

* 

P = P/Pf» 

* 2 

E = E/p^ 


We also need the following parameters, 
dynamics: 


Pr 


Re 


uRy 

k(Y - 1) 
u 


which occur extensively 


!!e 

k 


in fluid 


Here Pr and Re are abbreviations for the Prandtl number and Reynolds 
number respectively. Henceforth we drop the for notational 
convenience. For our numerical simulation we use Stokes' assumption: 


X = -2/3p. 


The nondimensionalized Navier-Stokes equations then take the form: 


Pj. + 


f \ . f ^ \ 49u 

(pu) + (pu + p)^ ^ 

3x 


(E^t 3Re 3x 3x^ (y - l)Pr Re 2 

oX 
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Rewriting the dissipative term (“D)q in terras of our new variable v 
we obtain: 


(C(v)v^)^ 



Our nondiraensionalized boundary conditions take the form: 


Supersonic Outflow 


8u 

9x 


0 


8x 


0 


Subsonic Outflow 


4 3u 
3Re 8x 





Subsonic Inflow 




8x 


0 


Y 3T 
Re Pr 8x 


T = where > 


(y - y + 2) 
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Supersonlc Inflow 


pu = 


4 3u 

3Re 3x " “2 


— ^ — — - a T = g 
Re Pr 3x 3 ®3 


where 


“2^2 


“3 > 


’1* 


We now propose a discrete version of these boundary conditions which are 
of extrapolation type and hence minimize the effect of the computational 
boundary. 


Supersonic Outflow 


^0 


T = T 
0 1 


We can connect qg to without any waves at all by putting 


PQ - Pl* 


Subsonic Outflow 

The boundary conditions proposed for supersonic outflow would work 
equally well for subsonic outflow. This corresponds to choosing = 0 in 
the subsonic outflow boundary conditions. If we choose ^2 ^ ^ 
discretized boundary conditions would be: 


4 

3Re 


~ 

Ax^ 


“2 = §1 
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T 


0 


= T 


1 


where we could put 


4 

3 Re 


(u^ - u^) 
Ax 


" “2 '' 2 * 


Then 


(4/3Re) - 02 Ax')u^ - Ax' 
^0 " (4/3Re) 


Clearly Uq Uj^ as Ax' 0. We choose by joining qj^ to qg by two 

wave interactions. 


Subsonic Inflow 


(pu)q = gi 


(pu>i - [(pu)2 - (pu) J 


u 


0 


u 


1 


Re Pr 


iVl 

Ax^ 




" “3 '^1 ^3 


(T 2 - 


Ti) 


Re Pr 


Ax 


“3 ^2* 


Here o^ > ( 7 ^ “ Y + 2)g^/2. 


Supersonic Inflow 

(pu)q = §1 = (pu)j - [(pu>2 - (pu) J 

4 S " ^0^ _ _ 4 ^^2 " 

IKE aF ~ “2 ''l ■ ®2 “ iRi- H " “2 "^2 

(T, - T.) (T - T,) 

Y 10 m___Y 21 

Re Pr aF “3 “ S 3 - Re Pr Ax ~ “3 4* 
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Here > (gj )/2, > g^. 

As mentioned earlier all the boundary conditions are such that ^ 
as Ax' 0 and they are a discretized version of the boundary conditions for 
the differential equations, i.e,, 

eR |a + Sq = g. 

It is easy to verify from (6,7) that for Ax' small enough we get bounded 
growth of the discrete version of the energy S(t) in time by choosing 
and appropriately. 

7 . NUMERICAL RESULTS 

Throughout the simulations we use the following parameter values: 

Pr = 0.7 
Y = 1.4 
Ax = 0 . 1 
Ax' = 0.000001 
At = 0.01. 

Simulation 1 : For the first simulation we use a very large value of the 

Reynolds number 

Re = 10^. 

We run the program with Rieraann initial data: 
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Pl = 1.0 

u = 3.0 

Xj 

Pl = 1.0 

and 

= 0.4734821 
u„ = 2.1393370 

K 

= 0.3333333. 

K 

At the left boundary our boundary conditions correspond to supersonic 
inflow. At the right boundary we have supersonic outflow boundary conditions. 

The solution to the Euler equations with initial data corresponding to 
the above Rieraann problem is a 3 shock moving to the right with a speed 
= 3.774567. 

The numerical simulations bear this out very well. The boundary 
conditions are radiative and allow the shock to pass through. Further for 
long time the solution stabilizes to a constant state. 
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RHO 



Figure 7.6 
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T=0.0 


u 



Figure 7.11 
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T=0.05 


U 



Figure 7.12 
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T=0.1 


U 



Figure 7.13 
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T=0.2 

U 



Figure 7.14 




Figure 7.15 
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Simulation 2 ; We use a low Reynolds number: 

Re = 500. 

Once more we run the program with Riemann initial data: 


Pl = 1.0 

\ 

Pl = 1.0 

and 

p_ = 1.625000 
R 

u„ = 0.3798263 
R 

p„ = 2.000000. 

At the left boundary we have subsonic inflow boundary conditions and at 
the right boundary the boundary conditions correspond to subsonic outflow. 

The solution to the Euler equations with the above Riemann initial data 
is a 1 shock moving to the left with a speed = 0,6124516, 

The numerical results, once again, have all the desirable properties we 
observed in Simulation 1. 
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T=0.0 


RHO 



Figure 7.16 
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T=0.8 


RHO 



Figure 7.19 
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T=0.0 


U 



Figure 7.21 
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T=0.2 


U 



Figure 7.22 






Figure 7.25 
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T=0.2 


P 



Figure 7.27 




Figure 7.28 




Figure 7*29 
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